/*********************************************************************
*
* Software License Agreement (BSD License)
*
*  Copyright (c) 2010, Heriot-Watt University, UK.
*  All rights reserved.
*
*  Redistribution and use in source and binary forms, with or without
*  modification, are permitted provided that the following conditions
*  are met:
*
*   * Redistributions of source code must retain the above copyright
*     notice, this list of conditions and the following disclaimer.
*   * Redistributions in binary form must reproduce the above
*     copyright notice, this list of conditions and the following
*     disclaimer in the documentation and/or other materials provided
*     with the distribution.
*   * Neither the name of the Heriot-Watt University nor the names of 
*     its contributors may be used to endorse or promote products 
*     derived from this software without specific prior written
*     permission.
*
*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
*  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
*  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
*  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
*  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
*  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
*  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
*  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
*  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
*  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
*  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
*  POSSIBILITY OF SUCH DAMAGE.
*
*  Author: Pedro Patron, Pierre-Yves Mignotte
*
*********************************************************************/

#ifndef _OSLCORE_VECTOR_H_
#define _OSLCORE_VECTOR_H_

#include <iostream>

#include <stdarg.h>
#include <stdint.h>

#include <osl_core/Iterator.h>

namespace osl_core {

typedef enum
{
	ONE, TWO, INF
} T_norm;

template<class T>
class Vector;

template<class T>
class Matrix;

template<class S>
Vector<S> operator +  (const Vector<S> & a, const S & k);

template<class S>
Vector<S> operator +  (const Vector<S> & a, const Vector<S> & b);

template<class S>
Vector<S> operator -  (const Vector<S> & a, const Vector<S> & b);

template<class S>
Vector<S> operator -  (const Vector<S> & a, const S & k);

template<class S>
Vector<S> operator -  (const Vector<S> & a);

template<class S>
Vector<S> operator *  (const S & k, const Vector<S> & a);

template<class S>
Vector<S> operator *  (const Vector<S> & a, const Vector<S> & b);

template<class S>
Vector<S> operator *  (const Vector<S> & A, const S & k);

template<class S>
S						dot (const Vector<S> & a, const Vector<S> & b);

template<class S>
Vector<S>			cross (const Vector<S> & a, const Vector<S> & b);

template<class S>
Vector<S> 			pot (const Vector<S> & a,   uint32_t n);

template<class S>
S						norm (const Vector<S> & a, T_norm norm = (T_norm) TWO);

/**
 * Allows change vector size dynamically
 * @class Vector
 */
template <class T>
class Vector
{	
	/**
	 * Negation
	 * \f$ x \in \Re^n \f$ 
	 * \f$ z = -x \to z_i = -x_i\f$
	 */
	template<class S>
	friend Vector<S> operator -  (const Vector<S> & x);
	
	/**
	 * Addition
	 * \f$ k \in \Re, a \in \Re^n \f$ 
	 * \f$ z = a + k \to z_i = a_i + k \f$
	 */
	template<class S>
	friend Vector<S> operator +  (const Vector<S> & x, const S & k);
	
	/**
	 * Addition
	 * \f$ a \in \Re, x \in \Re^n, y \in \Re^n \f$ 
	 * \f$ z = ax \to z_i = x_i + y_i \f$
	 */
	template<class S>
	friend Vector<S> operator +  (const Vector<S> & x, const Vector<S> & y);

	/**
	 * Substraction
	 * \f$ a \in \Re, x \in \Re^n, y \in \Re^n \f$ 
	 * \f$ z = ax \to z_i = x_i - y_i \f$
	 */
	template<class S>
	friend Vector<S> operator -  (const Vector<S> & x, const Vector<S> & y);

	/**
	 * Scalar Substraction
	 */
	template<class S>
	friend Vector<S> operator -  (const Vector<S> & a, const S & k);

	/**
	 * Scalar multiplication
	 * \f$ a \in \Re, x \in \Re^n \f$ 
	 * \f$ z = a*x \to z_i = a*x_i \f$
	 */
	template<class S>
	friend Vector<S> operator *  (const S & a, const Vector<S> & x);
	
	/**
	 * Vector multiplication of Hadamard product
	 * \f$ x \in \Re^n, y \in \Re^n \f$ 
	 * \f$ z = x.*y \to z_i = x_i*y_i \f$
	 */
	template<class S>
	friend Vector<S> operator *  (const Vector<S> & x, const Vector<S> & y);

	/**
	 * Cross product
	 * \f$ x \in \Re^n, y \in \Re^n \f$ 
	 */
	template<class S>
	friend Vector<S>	cross (const Vector<S> & x, const Vector<S> & y);
	
	/**
	 * Dot product 
	 * \f$ c \in \Re, x \in \Re^n, y \in \Re^n \f$ 
	 * \f$ c = x^T*y \to c = \sum_{i=1}^n(x_i*y_i) \f$
	 */
	template<class S>
	friend S		dot (const Vector<S> & x, const Vector<S> & y);

	/**
	 * Power 
	 * \f$ a \in \Re, x \in \Re^n \f$ 
	 * \f$ z = x^a \to z_i = x_i^a \f$
	 */
	template<class S>
	friend Vector<S> 	pow (const Vector<S> & x,   uint32_t a);

	/**
	 * Vector norm
	 * 
	 */
	template<class S>
	friend S				norm (const Vector<S> & a, T_norm norm );
	
	public:

		/**
		 * Constructor with initial value
		 * @param	numberElements number of components
		 * @param	initialValue	initial value
		 */
		Vector		(const uint32_t & numberElements = 3, const T & initialValue = (T) T() );

		Vector		(const uint32_t & numberElements, const uint32_t & maxElements, const T & initialValue);
		
		/**
		 * Copy constructor
		 */
		Vector		(const Vector<T> & source);

		/**
		 * Destructor
		 */
		virtual ~Vector ();
				
		/**
		 * Elements access
		 * @param	index index of the desired component
		 * @return reference to the component
		 */
		T &				 	operator []	(const uint32_t & index) const;

		/**
		 * Assignment
		 */
		Vector<T> &		operator =	(const Vector<T> & source);
		
		/**
		 * Vector size
		 * @note	ATTENTION!!! This function was previously incorrectly called length
		 */		
		uint32_t		size	() const;

		/**
		 * 
		 */
		uint32_t		allocate () const;

		/**
		 * Smallest element in the vector
		 */
		T 			minimum () const;

		/**
		 * Largest element in the vector
		 */
		T 			maximum () const;
		
		/**
		 * Vector length
		 * calculate from norm 2
		 */
		double	length	() const;

		/**
		 * Normalize vector
		 * @return vector with same origin and direction but unitary length
		 */
		Vector<T> 		normalize() const;
				
		/**
		 * Dynamic size modification
		 * @param numberOfElements	new dimension
		 *
		 * Modify vector to its new dimension
		 */
//		uint32_t		setSize				(const uint32_t & numberOfElements);

		/**
		 * Dynamic size modification with initial value
		 * @param numberOfElements	new dimension
		 * @param initialValue	initial value for new components
		 * Modify vector to its new dimension with new components with value initialValue
		 */
		 #if defined(WIN32) || defined(_WIN32)
		uint32_t		setSize				(const uint32_t & numberOfElements, const T initialValue = (T)0);
		#else
		uint32_t		setSize				(const uint32_t & numberOfElements, const T & initialValue = T());
		#endif
	
		
		uint32_t		setAllocate			(const uint32_t & numberOfElements);
		
		/**
		 * Set values to the received arguments
		 * The first parameters must match the vector size
		 */
		#if defined(WIN32) || defined(_WIN32)
		// Don't define this function under windows (until we fix the syntax)
		//void 		setValues (uint32_t n, ...);
		#else
		void 		setValues (uint32_t n, ...);
		#endif

		/**
		 * Scalar a x plus y
		 * \f$ a \in \Re, x \in \Re^n, y \in Re^n \f$ 
		 * \f$ y = a*x + y \to y_i = a*x_i + y_i \f$
		 */
		void  	saxpy (const T & a, const Vector<T> & x);
		
		/**
		 * Generalized scalar a x plus y
		 * \f$ a \in \Re, x \in \Re^n, y \in Re^n \f$ 
		 * \f$ y = A*x + y \to y_i = A(i,j)*x_i + y_i \f$
		 */
		void  	gaxpy (const Matrix<T> & a, const Vector<T> & x);
		
	protected:

		/**
		 * Buffer
		 */
		T *		m_data;

		/**
		 * Allocated size
		 */
		uint32_t	m_allocate;
		
		/**
		 * Dimension
		 */
		uint32_t	m_size;
};
			
/**
 * @param	v	vector
 * @param	i	index of first element
 * @param	j	index of second element
 * 
 * Interchange elements of index i and j
 */
template <class T>
void swap (Vector<T> & v, uint32_t i, uint32_t j);

/**
 * @param	v	vector
 *
 * Sort vector components by a direct selection method
 */
template <class T>
void sort (Vector<T> & v);

/**
 * Output vector stream
 * @param	out	output stream
 * @param	v	vector
 * @return	output stream with vector components
 */
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::ostream & operator << (std::ostream & out, Vector<T> & v);
#else
ostream & operator << (ostream & out, Vector<T> & v);
#endif


/**
 * Output vector stream
 * @param	out	output stream
 * @param	v	vector pointer
 * @return	output stream with vector components
 */
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::ostream & operator << (std::ostream & out, Vector<T> * v);
#else
ostream & operator << (ostream & out, Vector<T> * v);
#endif


/** 
 * Input vector stream
 * @param	in	input stream
 * @param	v	vector
 *
 * Modify vector to adjust to elements received from the input stream
 * The vector is adjusted to the number of received elements
 */
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::istream & operator >> (std::istream & in, Vector<T> & v);
#else
istream & operator >> (istream & in, Vector<T> & v);
#endif


/** 
 * Input vector stream
 * @param	in	input stream
 * @param	v	vector pointer
 *
 * Modify vector to adjust to elements received from the input stream
 * The vector is adjusted to the number of received elements
 */			
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::istream & operator >> (std::istream & in, Vector<T> * & v);
#else
istream & operator >> (istream & in, Vector<T> * & v);
#endif


/**
 * @param	v	vector as an ordered sequence of elements type T
 * @param	ele	element type T
 * @param	max maximum
 * @return index of the value 0<=k<=max that has v[0..k-1]<=ele<a[k]
 */
template <class T>
uint32_t binarySearch (const Vector<T> & v, T ele, uint32_t max);

/**
 * @param	v	vector
 *
 * Modify vector sorting its components by the heap method
 */
template <class T>
void heapSort (Vector<T> & v);

/**
 * @param	v	vector
 * @return	median value of the elements
 * 
 */
template <class T>
T median (const Vector<T> & v);


/**
 * @param	v	vector
 * @return the vector incremental values starting in value
 */
template <class T>
void incremental (const Vector<T> & v, const T & value = T());


/**
 * @param	v	vector
 * @return	number of unique elements in the vector
 */
template <class T>
uint32_t unique (const Vector<T> & v, Vector<T> * unique = (Vector<T>*)NULL, Vector<uint32_t> * index = (Vector<T>*)NULL);

/****************************************************************/
//						
//		VectorIterator
//						
/****************************************************************/

/**
 * Iterator protocol for vectors
 * @class	VectorIterator
 */
template <class T>
class VectorIterator : public Iterator<T>
{
	public:
			
		/**
		 * Constructor
		 */
		VectorIterator	(Vector<T> &);

		/**
		 * Copy constructor
		 */
		VectorIterator	(const VectorIterator<T> &);
		
		virtual bool	init		();
		virtual bool	operator !			();
		virtual bool	operator ++			();
		virtual bool	operator ++			(int);
		virtual T	operator ()			();
		virtual T &	operator =		 	(const T & newValue);
	
		/**
		 * Put the key in the previous current element
		 * @return true if there is current element
		 */
		bool		operator --	();	
		bool		operator --	(int);	

		/**
		 * Actual element index
		 */		 
		uint32_t	key		();
		
	protected:
	
		uint32_t	currentKey;
		Vector<T> &	theVector;
};


/****************************************************************/
//								
//		Vector Implementation
//								
/****************************************************************/

template <class T>
Vector<T>::Vector (const uint32_t & numberElements, const uint32_t & maxElements, const T & initialValue)
	: m_allocate(maxElements), m_size(numberElements)
{
	// initializes new vector
	assert((maxElements >= numberElements) && (numberElements > 0));
	
	// space allocation for the elements 
	m_data = new T[m_allocate];
	
	assert(m_data != 0);
	
	// initial value for each element
	for (uint32_t i = 0; i < m_size; i++)
		m_data[i] = initialValue;

}

//****************************************************************
template <class T>
Vector<T>::Vector (const uint32_t & numberElements, const T & initialValue)
	: m_allocate(numberElements), m_size(numberElements)
{
	// initializes new vector
	assert(numberElements > 0 );
	
	// space allocation for the elements 
	m_data = new T[m_allocate];
	
	assert(m_data != 0);
	
	// initial value for each element
	for (uint32_t i = 0; i < m_size; i++)
		m_data[i] = initialValue;
}


//****************************************************************
template <class T>
Vector<T>::Vector (const Vector<T> & source)
	: m_allocate(source.m_allocate), m_size(source.m_size)
{
	// copy constructor
	// space allocation for the elements
	m_data = new T[m_allocate];
	assert(m_data != 0);

	// copy source elements 
	for (uint32_t i = 0; i < m_size; i++)
		m_data[i] = source.m_data[i];

}

//****************************************************************
template <class T>
Vector<T> & Vector<T>::operator = (const Vector<T> & source)
{
	// asignation
	
	// free old space
	delete [] m_data;
	
	// create a new space for the elements
	m_data = new T[source.m_allocate];
	m_size = source.m_size;
	m_allocate = source.m_allocate;
	
	// copy source elements
	for (uint32_t i = 0; i < m_size; i++)
		m_data[i] = (T) source.m_data[i];
	
	return *this;
}

//****************************************************************
template <class T>
Vector<T>::~Vector ()
{
	// free dynamic allocated memory
	delete [] m_data;
	m_data = (T*)NULL;
	m_size = 0;
	m_allocate = 0;
}

//****************************************************************
template <class T>
T & Vector<T>::operator [] (const uint32_t & index) const
{
	// returns element of index index
	
	// checks index validity
	assert(index < m_size); 
	
	return m_data[index];
}

//****************************************************************
template <class T>
double Vector<T>::length () const
{
	// returns vector length

	return norm((*this), TWO);
}

//****************************************************************
template <class T>
uint32_t Vector<T>::size () const
{
	// returns vector size
	
	return m_size;
}

//****************************************************************
template <class T>
uint32_t Vector<T>::allocate () const
{
	// returns vector size
	
	return m_allocate;
}

//****************************************************************
template <class T>
T Vector<T>::minimum () const
{
	T minValue;

	if( m_size == 0 )
		minValue = 0;
	else
	{
		minValue = m_data[0];

		for (uint32_t i = 1; i < m_size; i++)
			if( m_data[i] < minValue ) 
				minValue = m_data[i];
	}
	return minValue;
}
//****************************************************************
template <class T>
T Vector<T>::maximum () const
{
	T maxValue;

	if( m_size == 0 )
		maxValue = 0;
	else
	{
		maxValue = m_data[0];

		for (uint32_t i = 1; i < m_size; i++)
			if( m_data[i] > maxValue ) 
				maxValue = m_data[i];
	}

	return maxValue;

}
//****************************************************************
template <class T>
Vector<T>  Vector<T>::normalize () const
{
	// returns normalized vector
	T n = norm((*this), TWO);

	if (n == (T) 0)
	{
		return (*this);
	}
	else
	{
		return (*this)*(1/n);
	}
}
//****************************************************************
template <class T>
uint32_t Vector<T>::setAllocate (const uint32_t & maxElements)
{
	assert(maxElements > 0);
	
	// create a new m_data area
	T * newData = new T[maxElements];
	assert(newData != 0);
	
	// copy all possible elements
	for (uint32_t i = 0; (i < maxElements) && (i < m_size); i++)
		newData[i] = m_data[i];

	// clean old m_data area
	delete [] m_data;

	// actualize m_data area	
	m_allocate = maxElements;
	m_data = newData;
		
	if (m_size > m_allocate)
	{
		m_size = maxElements;
	}

	// returns new dimension
	return m_allocate;
}

//****************************************************************
#if defined(WIN32) || defined(_WIN32)
//void Vector<T>::setValues (uint32_t n, ...)
#else
template <class T>
void Vector<T>::setValues (uint32_t n, ...)
{
	assert ((n == m_size) && ("Vector<T>::setValues() First argument must agree with vector size"));
	
	va_list ap;
	     
	va_start(ap, n);
	for (uint32_t i = 0; n > 0; i++, n--)
	{
		m_data[i] = va_arg(ap, double);
    }
	va_end(ap);
}
#endif

//****************************************************************
template <class T>
 #if defined(WIN32) || defined(_WIN32)
uint32_t Vector<T>::setSize				(const uint32_t & numberOfElements, const T initialValue)
#else
uint32_t Vector<T>::setSize				(const uint32_t & numberOfElements, const T & initialValue)
#endif
{
	// changes vector dimension
	
	assert(numberOfElements > 0);
	
	uint32_t oldSize(m_size);
	
	if (numberOfElements > m_allocate)
	{
		// create a new m_data area
		T * newData = new T[numberOfElements];
		assert(newData != 0);
	
		// copy all possible elements
		for (uint32_t i = 0; (i < numberOfElements) && (i < m_size); i++)
			newData[i] = m_data[i];

		// clean old m_data area
		delete [] m_data;

		// actualize m_data area	
		m_allocate = numberOfElements;
		m_data = newData;
	}
		
	m_size = numberOfElements;

	// if the size increases
	if (m_size > oldSize)
	{
		// set initial value to the new elements
		for (uint32_t i = oldSize; i < m_size; i++)
			m_data[i] = initialValue;
	}

	// returns new dimension
	return m_size;
}

//****************************************************************
template <class T>
void  Vector<T>::saxpy (const T & a, const Vector<T> & x)
{
	assert(m_size == x.size());
	
	for (uint32_t i = 0; i < m_size; i++)
	{
		m_data[i] = a*x[i] + m_data[i];
	}
}
		
//****************************************************************
template <class T>
void   Vector<T>::gaxpy (const Matrix<T> & A, const Vector<T> & x)
{
	assert(m_size == A.numberRows());
	assert(x.size == A.numberColumns());
	
	for (uint32_t i = 0; i < A.numberRows(); i++)
	{
		for (uint32_t j = 0; j < A.numberColumns(); j++)
		{
			m_data[i] = A[i][j]*x[j] + m_data[i];
		}
	}
}


/****************************************************************/
//							
//		VectorIterator Implementation
//
/****************************************************************/

//****************************************************************
template <class T>
VectorIterator<T>::VectorIterator (Vector<T> & aVector)
	: theVector(aVector)
{
		// iterator initialization
		init();
}

//****************************************************************
template <class T>
VectorIterator<T>::VectorIterator (const VectorIterator<T> & source)
	: theVector(source.theVector), currentKey(source.currentKey)
{
	// copy constructor
}

//****************************************************************
template <class T>
bool VectorIterator<T>::init ()
{
	// put actual key to the first element
	currentKey = 0;
	
	return operator ! ();
}

//****************************************************************
template <class T>
bool VectorIterator<T>::operator ! ()
{
	// returns true is there is current element
	return currentKey < theVector.size();
}

//****************************************************************
template <class T>
bool VectorIterator<T>::operator ++ ()
{
	// goes to the next elemente
	currentKey++;
	
	return operator ! ();
}

//****************************************************************
template <class T>
bool VectorIterator<T>::operator ++ (int)
{
	return ++(*this);
}

//****************************************************************
template <class T>
T VectorIterator<T>::operator () ()
{
	// returns current element
	return theVector[currentKey];
}

//****************************************************************
template <class T>
T & VectorIterator<T>::operator = (const T & newValue)
{
	// modify current element

	return theVector[currentKey] = newValue;
}

//****************************************************************
template <class T>
bool VectorIterator<T>::operator -- ()
{
	// back to the previous element
	if (currentKey != 0)
		currentKey--;
	
	return currentKey < theVector.size();
}

//****************************************************************
template <class T>
bool VectorIterator<T>::operator -- (int)
{
	return --(*this);
}

//****************************************************************
template <class T>
uint32_t VectorIterator<T>::key ()
{
	// returns current key
	return currentKey;
}

/****************************************************************/
//					
//	Vector functions
//					
/****************************************************************/

//****************************************************************
template <class T>
void swap (Vector<T> & v, uint32_t i, uint32_t j)
{
	T temp = v[i];
	v[i] = v[j];
	v[j] = temp;
}

//***************************************************************
template <class T>
void sort (Vector<T> & v)
{
	for (uint32_t i = v.size()-1; i > 0; i--)
	{
		// look for the position of the biggest element
		uint32_t largeposition = 0;

		for (uint32_t j = 1; j <= i; j++)
		{
			if (v[j] > v[largeposition])
			{
				largeposition = j;
			}
		}
		// place maximum in position i
		swap(v, i, largeposition);
	}
}


//***************************************************************
template<class T>
Vector<T> operator -  (const Vector<T> & x)
{
	Vector<T> z(x.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		z[i] = -x[i];
	}
	
	return z;
}

//***************************************************************
template<class T>
Vector<T> operator +  (const Vector<T> & x, const Vector<T> & y)
{
	assert(x.size() == y.size());
	
	Vector<T> z(x.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		z[i] = x[i] + y[i];
	}
	
	return z;
}
//***************************************************************
template<class T>
Vector<T> operator +  (const Vector<T> & x, const T & k)
{
	Vector<T> z(x.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		z[i] = x[i] + k;
	}
	
	return z;
}

//***************************************************************
template<class T>
Vector<T> operator -  (const Vector<T> & x, const Vector<T> & y)
{
	assert(x.size() == y.size());
	
	Vector<T> z(x.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		z[i] = x[i] - y[i];
	}
	
	return z;
}

//***************************************************************
template<class T>
Vector<T> operator -  (const Vector<T> & a, const T & k)
{
	Vector<T> z(a.length());
	
	for (uint32_t i = 0; i < z.length(); i++)
	{
		z[i] = a[i] - k;
	}
	
	return z;
}

//***************************************************************
template<class T>
Vector<T> operator *  (const T & a, const Vector<T> & x)
{
	Vector<T> z(x.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		z[i] = a*x[i];
	}
	
	return z;
}

//***************************************************************
template<class T>
Vector<T> operator *  (const Vector<T> & x, const T & a)
{
	return a*x;
}

//***************************************************************
template<class T>
Vector<T> operator *  (const Vector<T> & x, const Vector<T> & y)
{	
	assert(x.size() == y.size());
	
	Vector<T> z(x.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		z[i] = x[i] * y[i];
	}
	
	return z;
}


//***************************************************************
template<class T>
Vector<T>			cross (const Vector<T> & a, const Vector<T> & b)
{
	assert (a.size() == b.size());
	assert (a.size() == 3);
	
	Vector<T> s(3);

	s[0] = a[1]*b[2] - a[2]*b[1]; 
	s[1] = a[2]*b[0] - a[0]*b[2];
	s[2] = a[0]*b[1] - a[1]*b[0];

	return s;
}

//***************************************************************
template<class T>
T	dot (const Vector<T> & x, const Vector<T> & y)
{
	assert(x.size() == y.size());
	
	T c = T(0);
	
	for (uint32_t i = 0; i < x.size(); i++)
	{
		c = c + x[i]*y[i];
	}
	
	return c;
}

//***************************************************************
template<class T>
Vector<T> pow (const Vector<T> & x, uint32_t a)
{
	// optimised: not used of pow function

	Vector<T> z(x.size(),1);
	
	assert(x.size() == z.size());
	
	for (uint32_t i = 0; i < z.size(); i++)
	{
		for (uint32_t j = 1; j <= a; j++)
		{
			z[i] = z[i]*x[i];
		}
		// z[i] = pow(x[i], a);
	}
	
	return z;
}

//***************************************************************
template<class T>
T	norm (const Vector<T> & x, T_norm norm)
{
	T result = 0;
	
	switch (norm)
	{
		case ONE:
				for (uint32_t i = 0; i < x.size(); i++)
				{
					result = result + fabs(x[i]);
				}
			break;
					
		case TWO:
				result = sqrt(dot(x, x));
			break;
			
		case INF:
				for (uint32_t i = 0; i < x.size(); i++)
				{
					if (i == 0)
					{
						result = fabs(x[i]);
						continue;
					}
					if (result < fabs(x[i]))	
					{
						result = fabs(x[i]);
					}		
				}
			break;
			
		default:
			assert("Vector<T>::norm() Unknown norm called");
	}
	
	return result;
}

//****************************************************************
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::ostream & operator << (std::ostream & out, Vector<T> & v)
#else
ostream & operator << (ostream & out, Vector<T> & v)
#endif
{
	out << "( ";
	for (uint32_t i = 0; i < v.size(); i++)
		out << v[i] << ' ';
	out << ')';
//	out << " A:" << v.allocate();
//	out << " S:" << v.size();

	return out;
}

//****************************************************************
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::ostream & operator << (std::ostream & out, Vector<T> * v)
#else
ostream & operator << (ostream & out, Vector<T> * v)
#endif
{
	if (v)
		out << *v;
	else
		out << "()";

	return out;
}

//****************************************************************
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::istream & operator >> (std::istream & in, Vector<T> & v)
#else
istream & operator >> (istream & in, Vector<T> & v)
#endif
{
	uint32_t	i = 0;
//	uint32_t	len = v.size();
	char		ch;

	in >> ch;
	assert(ch == '(');
	
	in >> ch;
	
	while (ch != ')') 
	{
		in.putback(ch);
		in >> v[i];
		if (++i == v.size())
		{
			v.setSize(i + 10);
		}
		
		in >> ch;
	}

	// ajust input dimension
	if (i > 0)		
		v.setSize(i);

	return in;
}

//****************************************************************
template <class T>
#if defined(WIN32) || defined(_WIN32)
std::istream & operator >> (std::istream & in, Vector<T> * & v)
#else
istream & operator >> (istream & in, Vector<T> * & v)
#endif
{
	Vector<T> * p = new Vector<T>(100);

	in >> *p;
	v = p;

	return in;
}

//****************************************************************
template <class T>
uint32_t binarySearch (const Vector<T> & v, T ele, uint32_t max)
{	
	uint32_t low = 0;
	uint32_t high = max;

	while (low < high) 
	{
		// elements in v[0..low-1] are lower than ele
		// elements in v[high..max] are highers o equals than ele
 
		uint32_t mid = (low + high) / 2;

		if (v[mid] < ele)
			low = mid + 1;
		else
			high = mid;
	}

	// returns element that v[0..low-1] < ele <= v[low]

	return low;
}

//****************************************************************
template <class T>
void reorganize (Vector<T> & v, uint32_t low, uint32_t high)
{
	T	rootValue = v[low];
	uint32_t	position = low;

	while (position <= high) 
	{
		uint32_t childpos = 2 * position + 1;

		if (childpos <= high) 
		{
			if ((childpos + 1 <= high) && v[childpos+1] > v[childpos])
			{
				childpos++;
			}

			if (rootValue > v[childpos]) 
			{
				// we found the searched position 
				v[position] = rootValue;
				return;
			}

			// children to the actual position   
			v[position] = v[childpos];
			position = childpos;
		}
		else 
		{
			// position with no children 
			v[position] = rootValue;
			return;
		} 
	} 
}

//****************************************************************
template <class T>
void heapSort (Vector<T> & v)
{
	for (int i = v.size() / 2 - 1; i >= 0; --i)
	{
		reorganize(v, i, v.size() - 1);
	}

	for (uint32_t i = v.size() - 1; i > 0; --i)
	{
		swap(v, 0, i);
		reorganize(v, 0, i-1);
	}
}

//****************************************************************
template <class T>
T median (const Vector<T> & v)
{
	T  retValue;

	Vector<T> temp(v);

	for (int i = temp.size() / 2 - 1; i >= 0; --i)
	{
		reorganize(temp, i, temp.size() - 1);
	}

	for (uint32_t i = temp.size() - 1; i > 0; --i)
	{
		swap(temp, 0, i);
		reorganize(temp, 0, i-1);
	}

	if (temp.size()%2==0)
	{
		retValue = (T) (0.5*(temp[(uint32_t) (temp.size()/2 - 1)] + temp[(uint32_t) (temp.size()/2)]));
	}
	else
	{
		retValue = (T) temp[(uint32_t) (temp.size()/2)];
	}

	return retValue;
}

//****************************************************************
template <class T>
void incremental (const Vector<T> & v, const T & value)
{
	T	inc_val = value;
	for (uint32_t i = 0; i < v.size(); i++)
	{
		v[i] = inc_val++;
	}
}

//****************************************************************
template <class T>
uint32_t unique (const Vector<T> & v, Vector<T> * v_unique, Vector<uint32_t> * index)
{
	T		firstVal;
	uint32_t  retValue;
	T		rootValue;
	uint32_t	rootIndex;
	uint32_t	position;
	uint32_t	high;
	uint32_t	childpos;
	
	Vector<T> 		temp(v);
	Vector<uint32_t>	temp_index(temp.size());
	
	incremental(temp_index);
	
	// Sort the vector following the heap sort method
	//heapSort(temp);
	for (int i = temp.size() / 2 - 1; i >= 0; --i)
	{
		// reorganize(temp, i, temp.size() - 1);
		rootValue = temp[i];
		rootIndex = temp_index[i];
		position = i;
		high = temp.size()-1;
		
		while (position <= high) 
		{
			childpos = 2 * position + 1;
			if (childpos <= high) 
			{
				if ((childpos + 1 <= high) && temp[childpos+1] > temp[childpos])
				{
					childpos++;
				}
				if (rootValue > temp[childpos]) 
				{
					// we found the searched position 
					temp[position] = rootValue;
					temp_index[position] = rootIndex;
					break;
				}
				// children to the actual position   
				temp[position] = temp[childpos];
				temp_index[position] = temp_index[childpos];
				position = childpos;
			}
			else 
			{
				// position with no children 
				temp[position] = rootValue;
				temp_index[position] = rootIndex;
				break;
			} 
		} 
		
	}

	for (uint32_t i = temp.size() - 1; i > 0; --i)
	{
		swap(temp, 0, i);
		swap(temp_index, 0, i);
		
		// reorganize(v, 0, i-1);
		rootValue = temp[0];
		rootIndex = temp_index[0];
		position = 0;
		high = i-1;
		
		while (position <= high) 
		{
			childpos = 2 * position + 1;
			if (childpos <= high) 
			{
				if ((childpos + 1 <= high) && temp[childpos+1] > temp[childpos])
				{
					childpos++;
				}
				if (rootValue > temp[childpos]) 
				{
					// we found the searched position 
					temp[position] = rootValue;
					temp_index[position] = rootIndex;
					break;
				}
				// children to the actual position   
				temp[position] = temp[childpos];
				temp_index[position] = temp_index[childpos];
				position = childpos;
			}
			else 
			{
				// position with no children 
				temp[position] = rootValue;
				temp_index[position] = rootIndex;
				break;
			} 
		} 
	}
	
	// Count the number of values that are unique in the vector
	retValue = 0;
	for (uint32_t i = 0; i < temp.size(); )
	{	
		// Count unique values
		retValue++;

		// Jump duplicates
		for (firstVal = temp[i]; (i < temp.size() && temp[i] == firstVal); i++);
	}

	// If a vector is received to store the unique values...
	if (v_unique != NULL)
	{
		// Set the vector to the correspondent size
		v_unique->setSize(retValue);

		// If a vector is received to stored the index of those values
		if (index != NULL)
		{
			// Set the vector to the correspondent size
			index->setSize(retValue);
		}
		
		// Do the loop again but adding values to the vectors
		retValue = 0;
		for (uint32_t i = 0; i < temp.size();)
		{	
			// Count unique values
			(*v_unique)[retValue] = temp[i];
			if (index != NULL)
			{
				(*index)[retValue] = temp_index[i];
			}
			retValue++;
			
			// Jump duplicates
			for (firstVal = temp[i]; (i < temp.size() && temp[i] == firstVal); i++);
		}
	}
	
	return retValue;
}

} // namespace

#endif
